The Fractional Fourier Transform and Applications

NASA-CR-203439

NASA

National Aeronautics and
Space Administration

Ames Research Center

Moffett Field, California 94035

ARC275<R«vFab81)

7/v-&Y-c:<e

a/ & B4^

The Fractional Fourier Transform and Applications

David H. Bailey and Paul N. Swarztrauber

RNR Technical Report RNR-90-004

The Fractional Fourier Transform and Applications

David H. Bailey and Paul N. Swarztrauber

March 9, 1990

Abstract

The paper describes the "fractional Fourier transform", which admits computation by
an algorithm that has complexity proportional to the fast Fourier transform algorithm.
Whereas the discrete Fourier transform (DFT) is based on integral roots of unity e~ 2l "/ n ,
the fractional Fourier transform is based on fractional roots of unity e~ 2ir * a , where a is
arbitrary. The fractional Fourier transform and the corresponding fast algorithm are useful
for such applications as computing DFTs of sequences with prime lengths, computing
DFTs of sparse sequences, analyzing sequences with non-integer periodicities, performing
high-resolution trigonometric interpolation, and detecting signals with linearly drifting
frequencies. In many cases, the resulting algorithms are faster by arbitrarily large factors
than conventional techniques.

Bailey is with the Numerical Aerodynamic Simulation (NAS) Systems Division at
NASA Ames Research Center, MofFett Field, CA 94035. Swarztrauber is with the National
Center for Atmospheric Research, Boulder, CO 80307, which is sponsored by National Sci-
ence Foundation. This work was completed while Swarztrauber was visiting the Research
Institute for Advanced Computer Science (RIACS) at NASA Ames. Swarztrauber's work
was funded by the NAS Systems Division via Cooperative Agreement NCC 2-387 between
NASA and the Universities Space Research Association.

1

1. Introduction

The conventional fast Fourier transform (FFT) algorithm is widely used to compute dis-
crete Fourier transforms (DFTs) and discrete convolutions, and to perform trigonometric
interpolation. However, in some applications of the FFT, either the input is only partially
nonzero, or only part of the DFT result is required, or both. Nonetheless, the FFT algo-
rithm is ordinarily used unless the desired results can be more efficiently computed directly
from the definition of the DFT. We present here a technique that permits many of these
applications to be computed more efficiently. This same technique can also be applied in
other situations that do not admit efficient solution using standard FFTs.

The central concept here is a generalization of the DFT that is termed the fractional
Fourier transform (FRFT). It is defined on the m-long complex sequence x = (xj, < j <
m) as

m-l

G k (x,a) = J2^~ 2iriJhQ (1)

i=o

The parameter a will not be restricted to rational numbers and in fact may be any complex
number. Although this transform is defined for all integers /?, we will usually compute it
for the first m nonnegative values, i.e. < k < m. Note that the ordinary DFT and its
inverse are special cases of the fractional Fourier transform:

m-l

F k (x) = £ xje- 2 ^™ (2)

3=0

= G fc (aj, 1/m) < k < m (3)

-j m— 1

■**» = -£*;^ jVm (4)

= — G k (x, -1/m) < k < m (5)

TO v '

The discrete Laplace transform can also be written in terms of the fractional Fourier
transform.

If a is a rational number, the FRFT can be reduced to a DFT and can thus be evaluated
using conventional FFTs. Suppose that a = r/n, where the integers r and n are relatively
prime and where n > m. Let p be the integer such that pr = 1 (mod n). Then (1) can
be written

m-l

G k (x,a) = Y, x i*~ 27rijkrfn (6)

m— 1

= £ a ^ e " 2lri(pi)Wn (7)

3=0

m-l

^ E^e- 2,rWn (8)

3=0

= Fk(y) 0<k<n (9)

where subscripts are interpreted modulo n and where y is the n-long sequence (xo, «p> fyp,
••■, aj( m _i)p,0, 0, 0, •••,0). Thus the first n values of the FRFT can be computed by
performing an n-point FFT on the sequence y. We will take 5nlog 2 n as the cost of this
operation, since that is the number of floating point operations in a radix-2 FFT of size n.
We will see in section 5 that the cost of computing an ?i-point DFT of a sparse sequence
such as y can be reduced to about 5nlog 2 m floating point operations by employing a
decimation scheme. Also, we will see in section 3 that when n — qm, the first m values
of a DFT can be computed in roughly 5n log 2 m floating point operations by employing a
different decimation scheme. Either of these techniques would reduce the cost of computing
the FRFT using (9).

2. The Fast Fractional Fourier Transform Algorithm

An impetus for studying the fractional Fourier transform is the existence of an algorithm
for computing it that is significantly more efficient than the ones described in the previous
section. The computational cost of this algorithm is only about 20mlog 2 m floating point
operations, which is independent of the value of a. In particular, this cost does not depend
on whether or not a is rational or even real. The algorithm is based on a technique
originally due to Bluestein [6] and is related to what is known in the signal processing field
as the "chirp z-transform" (see [11] and [12]).

This algorithm can be derived by noting that 2jk = j 2 + k 2 — (k — j) 2 . The expression
for the FRFT then becomes

m-l

G h (*,a) = £ a . J . c -«& a +*-(*-i) a ]« (10)

i=o

m-l

= e

— irik 2 a V^ v.f,-^^ m{k-3) 2 a

Y, Xj-e"** "e***-" a (11)

m-l

= e-^-^yjz^ (12)

3=0

where the m-long sequences y and z are defined by

yj = x^ 2 " (13)

*j = e^' 2 " .(14)

Since the summation (12) is in the form of a discrete convolution, it suggests evaluation
using the well known DFT-based procedure. However, the usual DFT method evaluates
circular convolutions, wherein Zk-j = Zk-j+m' This condition is not satisfied here, but
instead Zk-j = Zj-k when k — j < 0. Fortunately, it is possible to convert this summation
into a form that is a circular convolution. Suppose we wish to compute the the FRFT for
< k < m. First select an integer p > m — 1, and extend the sequences y and z to length
2p as follows:

y 3 ; = m<j<2p (15)

Zj f = m<j<2p~m (16)

Zj = e «0'- 2 P) 3 2p-m<j<2p (17)

It can be seen by inspection that the first m values of Gk(x,a) satisfy

2p-l

G k (x,a) = e-**" £ yi z fc _,. < k < m (18)

It can also be seen by inspection that the sequence z now satisfies the required property
for a 2p-point circular convolution. Thus it follows that 2p-point DFTs may be used to
evaluate (18):

G k (x,a) = e-^^F^iw) 0<k<m (19)

where w is the 2p-long sequence defined by w^ — Fk(y)Fk(z). It should be emphasized
that this equality only holds for < k < m. The remaining 2p — m results of the final
inverse DFT are discarded. These three DFTs can of course be efficiently computed using
2p-point FFTs (for discussions of computing FFTs, see [1], [4], [5], [7], [9], [11], [16] and
[17]). ^ _ .

To compute a different m-long segment (?*+,(«, a), < k < m, it is necessary to
slightly modify the above convolution procedure. In this case z is as follows:

(20)
(21)
(22)

The remainder of the algorithm is unchanged. This complete procedure will be referred to
as the fast fractional Fourier transform algorithm.

The technique of converting the summation (12) into a circular convolution can also
be understood as the embedding of a Toeplitz matrix into a larger circulant matrix, which
admits evaluation using an FFT. Readers who wish to study the matrix formulation of this
algorithm are referred to [15].

We will assume that 5qlog 2 q floating point operations are required for a <j-point FFT.
This is the cost of the commonly used radix-2 complex FFT. If q is not a power of two^
the cost is somewhat higher, depending on the factors of q. If m is a power of two, the
obvious choice for p is p = m. Then the total computational cost of the above algorithm
is 20m log 2 m + 44m floating point operations, assuming that the exponential factors in
(19) and the FFT of the z sequence have been precomputed. Note that this cost is
approximately four times the cost of the usual m-point FFT.

The efficiency of the above algorithm depends on the efficiency of the underlying circular
convolution algorithm used to evaluate (18). A variety of fast convolution algorithms exist
in the literature (see [2] and [11]), and one of these may be more attractive than using
FFTs, depending on the computer hardware and the application. Further, FFTs over
fields other than the complex numbers (see [11]) may be profitably employed in some

z • = e ir *0'+*) 2 °

< j < m

zj =

m < j < 2p — m

z . — e ^(i+*-2p) 2 o(

2p — m < j < 2p

instances. However, for simplicity in assessing the computational costs of algorithms, we
will assume in the following that ordinary complex FFTs are employed. We will also assume
when assessing costs that the sizes of DFTs and FRFTs are powers of two, although the
algorithms to be discussed are valid for any size. Hereafter only the dominant term of these
computational costs will be listed, since lower-order terms are generally overshadowed by
implementation details, especially for large sizes.

3. Computing Longer and Shorter Segments of Results

Formula (19) gives the first m values of the fractional Fourier transform. Multiple
segments of m values can be computed even more efficiently. For example, suppose one
wishes to compute the first r results of an FRFT, where r = qm. These can be efficiently
computed by repeating formulas (19) through (22) for s = 0, m, 2m, •••,(<? — l)m. In
this case it is reasonable to assume that for each value of s, the DFT of the appropriate
z sequence has been precomputed. Significantly, the DFT of the y sequence needs to be
computed only once. Thus each segment of m results can can be computed by multiplying
the DFT of the y sequence by the precomputed DFT of the appropriate z sequence,
performing an inverse 2p-point FFT on the result, and multiplying by exponentials as in
(19). The dominant cost of computing r results in this manner is lOr log 2 m floating point
operations (when r is large compared to m).

The discussion in the previous paragraph assumed that one is computing the first r
values of the FRFT. However, any set of r results that is the union of m-long contiguous
segments can be computed in this manner, and the computational cost is the same. It is
not even necessary to assume that the starting indices s of the segments are multiples of
m.

Segments of the FRFT shorter than m in length can be efficiently computed by using
a decimation scheme similar to that used for computing a short segment of a DFT. Since
the fast algorithm for a short segment of a DFT is not widely known, it will be presented
first.

Suppose that one wishes to compute the first m values of the DFT of the n-long sequence
a?, where n = qm. We can write

m— 1 q~ 1

AM = EEw- w (' + ")"" (23)

1=0 j=0
g-1 m-1

= Y, e ~ 2irijk/n Y, x 3+i<i € ~ 2irilk/rn (24)

j=0 1=0

q-1

= J e -*^( yi ) 0<k<m (25)

j=o

where y^ denotes the m-long sequence (a?j+/g, < / < m). This means that the first m
values of the n-point DFT of x can be computed by writing the input as a q x m matrix
in column major order, performing m-point FFTs on each row, multiplying the results by

certain exponentials and then summing each column. The dominant computational cost of
this procedure is 5nlog 2 m floating point operations, which is less than that of computing
a full-sized n-point FFT whenever m < n/2.

Formula (25) may be easily generalized to compute m values of the DFT starting at
any index s. One then obtains

g-l

**+.(«) = X; e " 2 " i(fc+5)/n ^( z i) 0<k<m (26)

i=o

where Zj = (x j+ i q e~ 2iril ^ m , < I < m).

The scheme described by equation (25) can be seen as an abbreviated form of an FFT
algorithm sometimes called the "four step FFT" (see [1] p. 150, [4], [9] p. 569, and [17] p.
202 - 203). This FFT algorithm has received renewed interest lately due its suitability for
some advanced computer architectures, particularly systems with hierarchical memories.
Formula (25) is also closely related to the "discrete Zak transform" [3].

This same general method can be applied to compute an m-long segment of an n-point
FRFT:

m— lg— 1

©»(•,«) = £ 2>, +( , e - 2 * i « +, «> te (27)

1=0 j=0

= Ee"** ta I> i+ « i e- M,h " (28)

j-0 t=o

g-l

= Se-'^G^ga) < k < m (29)

i=o

where y^ = (xj + i q , < / < ra). As above, this means that the first m values of the
n-point FRFT of x can be computed by writing the input as a q x m matrix in column
major order, performing m-point FRFTs on each row, multiplying the results by certain
exponentials and then summing each column. Note that the exponential factors in (29)
may be incorporated into the exponential factors in formula (19) and thus do not need to
be counted in the computational cost. Hence the dominant cost of this scheme is 20nlog 2 m
floating operations, which is always less than that of directly computing an n-point FRFT-r
The corresponding formula for an m-long segment of an n-point FRFT beginning at
some index s is

g-l

G k+$ (x,a) = Y, e ~ 2 ^ Jik+ ' )a Gk{z jiq a) 0<k<m (30)

j=o

where Zj = (xj + i q e~ 2mlq ' a i < / < m). Note that the exponential factors in these Zj
vectors can be incorporated into the factors in (20) and (22) and thus do not need to be
counted. Hence the dominant cost of this scheme is 20nlog 2 m floating point operations,
the same as above.

4. Computing DFTs of Sequences with Prime Lengths

It is popularly believed that DFTs can only be evaluated by a fast algorithm when the
size m of the transform is a highly factorable integer, such as a power of two. But since
the DFT is merely a special case of the FRFT, it follows that the fast fractional Fourier
transform algorithm can be used to compute m-point DFTs even when m is a prime. The
computational cost is only about four times the cost of a power of two FFT of comparable
size. This algorithm was first presented in Bluestein's paper [6].

It should be noted that Bluestein's algorithm permits DFTs of any size (not just primes)
to be computed using power of two FFTs, which are generally the most efficient FFTs
available on computer system libraries. Thus depending on the computer system and
implementation, Bluestein's algorithm may be more efficient for computing DFTs of certain
sizes than conventional FFTs. On some computer systems, such as hypercubes, Bluestein's
algorithm is preferable for computing DFTs of any size that is not a power of two [15].

Perhaps one reason that Bluestein's algorithm has not received more attention is that
another FFT algorithm for prime m was discovered at about the same time by Rader [13].
This algorithm reduces a DFT of prime size m to a circular convolution of size m — 1,
which in theory is only about half as costly as Bluestein's algorithm. Unlike Bluestein's
algorithm, however, Rader's algorithm does not generalize to the case of arbitrary a, and
there is little flexibility in the permissible sizes of FFTs that may be used for efficient
evaluation.

There are a number of other important applications of the basic technique behind
Bluestein's algorithm, which technique we have termed the FRFT. Several of these appli-
cations will be presented here, including some that were originally presented in [12],

5. Computing DFTs of Sparse Sequences

In this section we will describe efficient algorithms for computing DFTs of sequences
with large numbers of zero elements. DFTs of such sequences arise in trigonometric inter-
polation, as we shall see in section 7. Another common application is the computation of
the "linear" or "aperiodic" convolution, which is defined on the n-long sequences x and y
as the 2n-long sequence z — (£j=o ZjVk-jy < & < 2n — 1), where the subscript k — j
is restricted to the range < k — j < n — 1. Linear convolutions arise in applications
as diverse as multiprecision arithmetic and the numerical solution of partial differential
equations. This type of convolution is normally evaluated by padding the n-long input
sequences as and y with zeroes to length 2n, and then computing their 2n-point circular
convolution using FFTs.

Three cases of sparse DFTs will be considered. First we will compute the complete
DFT of an input sequence where only an initial segment is possibly nonzero. Second, we
will compute segments of such a DFT. Third, we will compute segments of the DFT of an
input sequence that has only a few nonzero segments.

Let n = qm. Suppose we wish to compute the DFT of an n-long sequence as that
is known to be zero except for the first m entries. This can be done by employing a

conventional FFT- based decimation scheme, as follows:

m-l

W») = E x,e- 2 "«* + '')/" (31)

E **•«"

i=o

-2*ij(k+lq)/n

m-l

E *i c '

-2irijk/n e -

-2-Kijlfm

fl(»»)

0</

< 771

(32)

< ft < q (33)

where #*. = (a^e -2 **'*/", < j < m). Thus for each &, < k < q, the m values of the
n-point DFT beginning at index k and increasing with stride q can be computed with a
conventional m-point FFT on the input sequence multiplied by certain exponentials. The
dominant cost of this procedure is 5nlog 2 m floating point operations, which is less than
that of an n-point FFT whenever m < n/2. This scheme will be used in section 7.

If only a single contiguous segment of the n-point DFT is required, then this approach
is unsatisfactory, since one can only obtain DFT results with a spacing of q from the above.
Further, the conventional method for computing a segment of a DFT in section 3 is also
unsatisfactory since each of the DFTs in the summation (26) must still be evaluated, with
no savings in computation. However, the FRFT can be effectively applied to this situation.

The simplest case is when only m values of the n-point DFT are required. In this case
we have

m-l

F k+t (x) = £ ^c-*****)/- (34)

3=0

= Gfc+,(y,l/n) 0<k<m (35)

where y = (xj, < j < m). The dominant cost here is the cost of an m-point FRFT,
or 20m log 2 m operations. This is less than the cost of conventional methods whenever
m < n/4 (for sufficiently large n).

Multiple sections can be computed even more economically by merely repeating (35)
for multiple values of a. As indicated in the first paragraph of section 3, savings result
because the FFT of the sequence y in (18) needs to be computed only once. The dominant
cost of computing r elements in this manner is lOr log 2 m floating point operations (when
r is large compared to m). Thus up to ra/2 elements of the DFT may be computed by
this method, and the cost is still less than computing all or part of the n-point DFT by a
conventional method (for large n).

We will now examine the case in which more than one m-long segment of the input
sequence x is possibly nonzero. An m-long segment of the n-point DFT of x can be written

F h+ ,(x) = £2>, + , m e- 2 " 0+ ' m)( * + * )/n (36)

j=o i=o

q-i m-l

= E € _2 ^ fe+ ^ m/n ]T x j+ i m e- 2 * ij ( k+a)/n (37)

f=0 j=0

8

9-1

= H e ~ 2wii(k+ ' )l/q Gk + .(ViA/n) 0<k<m (38)

1=0

where y t = (aj J+ / m , < j < m). With the assumption that only a few m-long input
segments are zero, it follows that most of the terms of this summation are zero. Thus only
a few of these FRFTs need to be actually computed.

We can assume that the exponential factors in (38) have been incorporated into the
exponential factors of formula (19). Suppose that only t of the-input segments y t contain
possibly nonzero elements. Then the total computational cost of computing an m-long
segment of the DFT is on the order of 20£mlog 2 m operations. This is less than the cost
of computing a segment of a DFT using the conventional scheme (33) whenever t < <jr/4
(for large n).

As indicated in section 3, multiple segments of results can be even more efficiently
computed by repeating (38) for several values of s. The dominant cost of computing r
values of the DFT in this way is lOtr log 2 m floating point operations (when r is large
compared to m). This is less than conventional methods whenever tr < n/2 (for large n).
Thus this method is not appropriate unless only a few segments of the input are possibly
nonzero and only a few segments of the result are desired.

6. Analyzing Sequences with Non-Integer Periodic Components

In many practical applications of DFTs and FFTs, notably in signal processing, it is
not really true that the underlying time series being studied is exactly periodic with period
m, where m is the size of the data sequence. Instead, the sampling interval and sample
size are usually selected arbitrarily for the convenience of hardware and software. Thus
in such applications it is the rule rather than the exception that a data sequence contains
components with periods that are not exactly whole numbers.

One instance of such an application is in computing solar positions. Daily solar altitude
and declination angles are readily available from almanacs. Such positions are not periodic
with period 365 days, but they are very nearly periodic with a fractional period of approx-
imately 365.2422 days. In other words, the apparent motion of the sun has a fairly simple
Fourier representation based on a period of 365.2422 days, but the representation is not
simple based on a period of 365 days. In this example we know from a priori information
what the true period is. But it general we usually do not even know that. Thus we see£
techniques that can determine such unknown periods and compute spectra based on the
true period once it is known.

Let x be the m-long sequence (e 2 **^ m , < j < m) where /3 is unknown and not an
integer. Then the DFT of x is

m-l
**(») = !£ e 2^i/3/m e -2iriiVm (39)

m-l

= Yl e 2irij ^ -/e)/m (40)

1 _ e 2irtj3

= , . (a ... 0<k<m (41)

Thus the spectrum of such a sequence is not a "spike" at a single index but instead is spread
over all indices. The magnitude of Fjt(x) is greatest for k closest to /3 and declines rapidly
with increasing distance from {3. It is fairly easy to determine the unknown frequency /3 to
within one or two units in this manner, but it is difficult to determine its true value more
accurately with such conventional techniques.

The fractional Fourier transform is well suited to studying these sequences because the
a in its definition may be arbitrary. In particular, one can use the FRFT to efficiently
compute a high resolution spectrum near /3. Let b be the greatest integer less than /?. Then
m spectra values starting at index b and increasing with interval S may be computed as
follows:

m-l

Fb+u(m) = E ^c-*^*")/** (42)

i=o

m-l
= E a iC -a*fcVm c -2iryM/m ^

= G h (y,6/m) 0<k<m (44)

where y = (xje~ 2m ^ m y < j < to). If fewer than to or greater than to spectral values
are desired, these can be efficiently computed by applying the techniques of section 3.

Suppose now that k is the index of the largest magnitude of this high resolution spectra.
Then a more accurate estimate of the true period /? is 6 + kS. With this new estimate
of /3, one can proceed to compute a period-based spectra. For example one could then
set 7 = (b + k6)/b and then apply formula (44) to compute the spectrum at the points
0, 7, 27, • • • , (m — 1)7. This spectrum will now have a fairly sharp peak at location 67.
This two-step procedure may be repeated to produce even more accurate estimates of j3
and even sharper spectra based on these estimated periods.

7. Performing High-Resolution Trigonometric Interpolation

The FFT is frequently used to perform trigonometric interpolation. For example, given
a sequence x in tabular form, we may wish to determine approximate midpoint values -
Xj+i/2* In the previous sections we have used the "aliased" form of the DFT:

m-l

F k (x) = X> ie - 2irWm 0<k<m (45)

i=o

In this section we will use the "unaliased" form of the DFT, which is defined on x =
(xj, — to/2 < j < m/2) as

m/2-1
A(«) = E Xje- 2 ™* 1 ™ -m/2<k<m/2 (46)

3= -m/ 2

10

The inverse unaliased DFT is defined analogously. For most purposes, the unaliased form
is not as convenient as the aliased form, because two different formulas are required for
m even and m odd (we will assume in the following that m is even). Further, most FFT
software for computing DFTs assumes the aliased form.

In discussions of trigonometric interpolation, however, the unaliased form is preferred.
The problem is that although Fk(x) = Fk{x) for each fe, their interpolated values may
be different (depending on how they are defined), so that Ffe +1 / 2 («) ^ i ? fc +1 / 2 (aj). This
disagreement stems- from the fact that the function /*(£) = e %kt has an infinite number
of alternate characterizations (or "aliases") on the points tj — 2nj/m, namely /fc+ m (£) =
e t(k+m)t £ or gjj m tegers m. These functions are indistinguishable from one another on the
points tj. However, at intermediate or interpolated values fk(t) and fk+ m (t) are different.
Therefore a choice of basis functions for the trigonometric expansion of a function /(£)
must be made. The most common choice is the set (/*(£) — e***, —m/2 < k < ra/2). This
set has minimum variance on < t < 2ir and thus ensures that the interpolated values will
be as smooth a possible.

For computational purposes, the unaliased form can be easily reduced to the aliased
form. One approach is to define m-long sequences y and Y such that y$ — (— l) J 'x,-_ m /2
for < j < m and Y k = (-1)* F k _ m/2 (x) for < k < m. Then

Y k = (-l) k F k -m,2(*) (47)

m/2-l

= e"^* J] x ie _2lrii(fc ~ m/2)/m (48)

j=— m/2
m-1

= E *;-m/ 2 e- 2 ^- m/2 > fc/m (49)

i=o

m-l

- E^ e " 2 " iVm (50)

= Fh(v) 0<k<m (51)

Thus Y is the ordinary m-point DFT of y. Alternately, one can define yj = Xj for
< j < m/2 and yj = a:j_ m for m/2 < j < m, and similarly define Yk in terms of Fk.
Then again Y is the ordinary DFT of y. Either way, FFT software designed to computer
the usual m-point aliased DFTs can be readily employed to compute m-point unaliased
DFTs.

We will now develop the traditional approach to trigonometric interpolation. Suppose
we wish to compute q interpolated values between each of the m values of x. Let n = qm.
We first compute the inverse unaliased DFT of x and then extend this spectra to the n-long
sequence Y as follows:

Y k = F k -\x) \k\ < m/2 (52)

Y- m/2 = \P: m/2 (x) (53)

11

m

/2 = 5^i /a (4) (54)

Y k = -n/2<k< -m/2 (55)

I* = m/2 <k<n/2 (56)

The splitting of ■FT m / 2 (£) i n the above is necessary to produce the most accurate interpo-
lated values. For example, if this is not done then trigonometric interpolation on purely
. real data will yield complex values. This split is not necessary if m is odd.

The 7i-long vector y of interpolated results yj = Xj/ qy where n = qm, can now be
computed from the unaliased n-point DFT of Y:

m/2

Vi = E r fc e- 2 ™ (j/ * )A/m (57)

fc=-m/2
n/2-1

= E Y k e-**M n -n/2<k<n/2 (58)

h--n/2

This approach requires an m-point inverse FFT followed by an n-point forward FFT, so
that the dominant cost is 5n log 2 n floating point operations.

This cost can be reduced somewhat by applying a decimation technique similar to that
employed in equations (31) through (33):

m/2

WWt = £ ««-**<*+*)/« (59)

k——m/2

m/2
_ £ y fcC -2**,yn e -2«MB/n (gQ)

*=-m/2

m/2-1
— Y^ V, e -2irikj/n e -2triklq/n _|_ y -irimj/n e -*imlq/n /g^\

k=— m/2

- F / (Z i ) + (-l) / y m/2 e-^* (62)

= l|(W r i ) + (-iy2y m / 2 cos(7rj7 ? ) -m/2</<m/2 0<j<? (63)

where Z,- = (Y k e- 2rik ^ n , -m/2 < & < m/2) and where Wj = Zj except that (Wj)_ m / 2 ,
the leftmost element of Wj, is zero. The formulas (62) and (63) differ only in that the
second is more appropriate for working with real data, since a real-valued FFT may be
used for computation. The dominant cost of computing all n interpolated values in this
fashion is 5nlog 2 m floating point operations.

In many applications of trigonometric interpolation, only a segment of the n interpo-
lated values are required. In this case the cost can be significantly reduced by applying
the FRFT, using techniques similar to those in section 5. Suppose we wish to compute
only the first m interpolated values in (58). This expression can be converted to a FRFT

12

as follows:

» = En-^e-^-"^" (64)

k=0

m

= e"'^J2Yk- m /2e- 2 " k ' /n (65)

fe=0

= e^^(Z,l/n) 0<j<m (66)

where Z is the (m + l)-long sequence (lfe_ m / 2 , < k < m).

Suppose for a minute that m is a power of two, and that the most efficient FFT sizes
employed for evaluating the FRFT on a given computer are also powers of two. It might
at first appear that increasing the size of an FRFT from m to m + 1 would require that
the size of FFTs used in computing the FRFT must be doubled to the next highest power
of two, with a corresponding sharp increase in computational cost. However, recall that
the parameter p of (18) can be any integer greater than one less than the size of the
FRFT. Thus 2m-point FFTs can still be used to evaluate FRFTs of size m + 1, and the
computational cost is virtually unchanged.

More generally, we can compute an m-long segment z of interpolated values beginning
with zq = x, and increasing with interval 8 as follows:

m/2

*j = E Yke-^'+M™ (67)

k—— m/2
m

= E *i_ /ie -* - (»— ^X»*-»/» (68)

m

= e^' + ^^^-W2 e " 2,n ' Wme " 2,riiW/m (69)

k=o

= e^'^G^Z.S/m) 0<j<m (70)

where Z is the (m + l)-long sequence (Yk- m /2e~ 2irik ^ m , < k < m).

Since the exponential factors in (70) can be incorporated into the precomputed expo-
nentials in (19), we do not need to count the cost of these multiplications. The dominant
cost of the algorithm is thus only 25mlog 2 m floating point operations. This is less than,
the cost of the conventional procedure in equations (62) and (63) whenever m < n/A (for
large n). Larger sections of results can be efficiently obtained by repeating (70) for different
starting indices s as indicated in section 3.

8. Detecting Signals with Linearly Drifting Frequencies

The authors first learned of the problem of detecting signals with linearly drifting
frequencies in discussions with scientists associated with the Search for Extraterrestrial
Intelligence (SETI) project at NASA Ames. Subsequently, the authors have learned that
similar problems arise in such fields as high- accuracy position determination using the
global positioning system satellites. In fact, signals with drifting frequencies will arise in

13

any situation where there is acceleration in the direction between the emitter of the signal
and the receiver. It often can be assumed that over the period of observation the frequency
drift is linear, so we will make this assumption in the following discussion.

Suppose that y is a time series of length n = qm that is the sum of a signal with a
linearly drifting frequency plus Gaussian noise. We will assume that the frequency is less
than m. Thus the time series can be written

Vj = ae 2 ^ w+ir ^ m + 6^ 0<i<n (71)

where a and b are scaling factors, and g$ are complex Gaussian deviates with unit variance.
We will be interested in applications where is a is much less than 6, or in other words where
the amplitude of the signal is much less than the amplitude of the noise.

The usual approach to detecting and analyzing such signals is to consider the n-long
time series asagxm array in row major order. We first compute the rn-point DFT of
each row in this array (i.e. the DFT of successive contiguous m-long segments of y), and
then compute squared magnitudes of each resulting array entry (see [8] and [10]):

X jtk = \F k (zj)\ 2 0<j<q, 0<k<m (72)

where Zj = (yj+i m , < / < m). When the frequency of the signal is constant (i.e. when
r is zero), detection of the signal can be achieved by merely summing the columns of
X to form an m-long vector. The "bin" containing the signal will then be evident as a
statistically large entry in this sum vector.

When the frequency of the signal to be detected is not constant (i.e. when r is nonzero),
this method fails since the "bin" containing the signal is different from segment to segment,
and thus no entry of the sum vector will contain a statistically large value. Nonetheless,
the signal is evident in this array as a "line" with a certain slope, and thus in theory it is
possible to detect such signals by computing sums in the X array over lines with a range
of slopes a:

9-1

S*( a ) = E X i,*int(M-«j) 0<fe<m (73)

3=0

where nint denotes nearest integer. When k is close to w and a is close to 2mr, then 5fc(a)
has a statistically large value and detection occurs. Unfortunately, the direct computation
of all these sums over a range of possible slopes a is very expensive [8].

Consider what appears at first to be an even more expensive calculation, namely com-
puting sums of values interpolated along the rows of X. In particular, consider the array
of sums

S r , k = Y, x i*H«+*r)i 0<r<q, 0<k<m (74)

where the subscript k + (a + Sr)j is interpreted literally as its fractional value, namely
as the result of trigonometric interpolation along the rows of the X array. Except for the

14

fractional subscripts, this corresponds to computing the sums 5fc(a r ), where a, = a + St.
The location (r, k) of the maximum value in the 5 array gives the drift rate and the initial
frequency of the detected signal.

This array of sums can be efficiently computed using the FRFT as follows. Define Xj t k
for negative k so that Xjj, = X^k+m* Let (xjj, —m/2 <l< m/2) be the unaliased inverse
m-point DFT of the j-th row of X. Then we can write

m/2-1

*** = E »i,ie- 2 " W/m ■ (75)

I=-m/2
m/2-1

XsMl*** = E «««-*-«*+<-+*»•/- (76)

f=-m/2

Thus

St* = E E 3 W e- 2 *[* + < a+ *W> m (77)

j=0 f=-m/2

",-1

, — 2iriajl/m — 2iri6rjlfm

E^/e" ' "e

c -a«M/m ( 7g J

m/2-1

I=-m/2 Li=°
m/2-1

£ v-""' (79)

i=-m/2

where

s r>l = Gi(z h Sr/m) (80)

*i = {*ue-*** /m , < i < <?) (81)

Hence this summation algorithm consists of performing inverse FFTs on rows of the X
array, followed by FRFTs on the resulting columns (with certain exponential multipliers),
followed by forward FFTs on the resulting rows. The computational cost of entire detection
algorithm is approximately 10nlog 2 n floating point operations, or in other words only
about twice the cost of an n-point complex FFT. For complete details, see [14].

9. Conclusions

We have described a number of applications of the FRFT and its associated fast algo-
rithm. The flexibility and efficiency of these methods permits new solutions to a variety
of problems that involve discrete Fourier analysis. In many cases, the resulting algorithms
are faster than the conventional methods of solution that utilize ordinary FFTs. In fact,
for sufficiently large problem sizes, these methods are often faster than the conventional
algorithms by arbitrarily large factors.

15

Acknowledgment

The authors wish to thank Izidor Gertner of the CUNY Graduate Center for informative
discussions on this subject.

References

1. Agarwal, R. C, and Cooley, J. W., "Fourier- Transform and Convolution Subroutines
for the IBM 3090 Vector Facility", IBM Journal of Research and Development, vol.
30 (1986), p. 145 - 162.

2. Agarwal, R. C, and Cooley, J. W., "New Algorithms for Digital Convolution", IEEE
Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-25, no. 5
(1977), p. 392 - 410.

3. Auslander, L., Gertner, I. and Tolimieri, R., "The Discrete Zak Transform: Applica-
tions to Time- Frequency Analysis and Synthesis of Non-Stationary Signals", IEEE
Transactions on Acoustics, Speech and Signal Processing, to appear.

4. Bailey, D. H., "FFTs in External or Hierarchical Memory", Journal of Supercomput-
ing, vol. 4 (1990), p. 23 - 35.

5. Bailey, D. H., "A High- Performance FFT Algorithm for Vector Supercomputers",
Internationa] Journal of Supercomputer Applications vol. 2 (1988), p. 82 - 87.

6. Bluestein, L. I., "A Linear Filtering Approach to the Computation of the Discrete
Fourier Transform", IEEE Transactions on Audio and Elect roacoustics, vol. AU-18,
no. 4 (1970), p. 451 - 455.

7. Cooley, J. W., and Tukey, J. W., "An Algorithm for Machine Computation of Com-
plex Fourier Series", Mathematics of Computation, vol. 19 (1965), p. 297 - 301.

8. Cullers, D. K., Linscott, I. R., and Oliver, B. M., "Signal Processing in SETI",
Communications of the ACM, vol. 28 (1985), p. 1151 - 1163.

9. Gentleman, W. M., and Sande, G., "Fast Fourier Transforms - For Fun and Profit",
AFIPS Proceedings, vol. 29 (1966), p. 563 - 578.

10. Gertner, I., "Optimal Detection and Separation of Chirp Signals", Proceedings of
ICASSP-90, 1990, to appear.

11. Nussbaumer, H. J., Fast Fourier Transform and Convolution Algorithms, Springer-
Verlag, New York, 1981.

12. Rabiner, L. R., Schafer, R. W., and Rader, C. M., "The Chirp z-Transform Algorithm
and Its Application", Bell System Technical Journal, vol. 48 (1969), p. 1249 - 1292.

16

13. Rader, C. M., "Discrete Fourier Transforms When the Number of Data Samples Is
Prime", Proceedings of the IEEE, vol. 56 (1968), p. 1107 - 1108.

14. Swarztrauber, P. N., and Bailey, D. H., "Efficient Detection a Continuous Wave
Signal with a Linear Frequency Drift", submitted for publication.

15. Swarztrauber, P. N., Sweet, R. A., Briggs, W. L., Henson, V. E., and Otto, J.,
"Bluestein's FFT for Arbitrary N on the Hypercube", submitted for publication.

16. Swarztrauber, P. N., "FFT Algorithms for Vector Computers", Parallel Computing,
1 (1984), p. 45 - 63.

17. Swarztrauber, P. N., "Multiprocessor FFTs", Parallel Computing, vol. 5 (1987), p.
197 - 210.

17

I